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A new variant of the pencil-beam (PB) algorithm for dose distribution calculation for radiother- 
apy with protons and heavier ions, the grid-dose spreading (GDS) algorithm, is proposed. The GDS 
algorithm is intrinsically faster than conventional PB algorithms due to approximations in convo- 
lution integral, where physical calculations are decoupled from simple grid-to-grid energy transfer. 
It was effortlessly implemented to a carbon-ion radiotherapy treatment planning system to enable 
realistic beam blurring in the field, which was absent with the broad-beam (BB) algorithm. For a 
typical prostate treatment, the slowing factor of the GDS algorithm relative to the BB algorithm 
was 1.4, which is a great improvement over the conventional PB algorithms with a typical slow- 
ing factor of several tens. The GDS algorithm is mathematically equivalent to the PB algorithm 
for horizontal and vertical coplanar beams commonly used in carbon-ion radiotherapy while dose 
deformation within the size of the pristine spread occurs for angled beams, which was within 3 
mm for a single proton pencil beam of 30° incidence, and needs to be assessed against the clinical 
requirements and tolerances in practical situations. 
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I. INTRODUCTION 

Ceaseless efforts for accuracy improvement and con- 
stant progress in computing technology have made 
the pencil-beam (PB) algorithm be the standard 
method for dose distribution calculation in heavy 
charged particle radiotherapy with protons and heav- 
ier ionSiiiMi££&LMiI2iII In the PB algorithm, a treat- 
ment beam is divided into elementary pencil beams 
with developing transverse spread as they penetrate 
through heterogeneous medium to handle spatial mod- 
ulation of beam scatter that is ignored in the broad- 
beam (BB) algorithmic The dose distribution will be 
formed with superposition of the pencil beams using 
kernel-convolution techniques with variations in algorith- 
mic implementation, for example, in choice of the coordi- 
nate system, order of the multiple integrals, and numeri- 
cal approximations, which greatly influence the accuracy, 
speed, complexity, and generality of the code i 12 i 13 

Though the PB algorithm may be sufficiently accurate 
and fast for dose calculation in treatment planning in 
the present form, demand for faster calculation methods 
may always remain, for example, for optimization in the 
intensity-modulated radiotherapy with scanned charged 
particle beams,— ^ and for adaptive radiotherapy under 
image guidance ) 15 i 16 which will ultimately accommodate 
on-site re-planning for an immobilized patient quickly be- 
tween imaging and treatment. Pursuing faster computa- 
tional algorithms might be critical for the innovation to 
happen. 

This paper presents one of such approaches, where we 
briefly review the BB and PB algorithms, describe the 
new algorithm, demonstrate the effectiveness in carbon- 



ion radiotherapy, evaluate the accuracy with a modeled 
proton pencil beam, and discuss the applicability for 
scanned beams. 



II. MATERIALS AND METHODS 

A. The broad-beam algorithm 

In the BB algorithm, dose D at point f is resolved 
into the BB dose and the penumbra effect*^ The BB 
dose Dbb, or equivalently dose per fluence of the incident 
beam, is given either theoretically or experimentally as 
a function of water- equivalent depth w that is calculated 
with the ray-tracing integral of effective density p from 
the beam source r$ in radial direction v = (f— f\j)/\r — 
r*o | . The penumbra effect gradates the field edge with the 
error function erf (a:) = (2/y / 7r) J* e~ u2 du of the signed 
closest distance to the geometrical field edge, t (t > for 
r in the field, t < otherwise) , 



D(f) = D BB (w(7))~ 
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where the projected transverse spread Ot(r) is given ei- 
ther experimentally or theoretically. 

The formulation of the penumbra effect simulates the 
dose-collecting process at point f from uniformly and 
continuously distributed invariant Gaussian sub-beams 
in the field. The assumed uniformity and invariance of 
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the sub-beams restrict the validity of the model to a nar- 
row penumbra region where the local-homogeneity ap- 
proximation may be valid. 



B. The pencil-beam algorithm 

In convolution algorithms, a dose distribution is gen- 
erally calculated by kernel integral, 



D(r) = J T(p)h(p,r)d 3 p, 



(3) 



where T(p) is the total energy released per mass (terma) 
from the radiation at point p and kernel function h(p,r) 
is the terma fraction transfered to point In the PB 
algorithm, the terma equals the BB dose in the beam 
field or zero otherwise and is transversely spread by a 
planar Gaussian kernel, 
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^S({f-p)-v(p)), (5) 



where the Dirac S function restricts the spreading in the 
plane perpendicular to the PB direction v, leading the 
convolution to an areal integral of a pencil kernel in the 
fieldpii^ In the numerical integration, several tens or more 
terma-emitting points are usually arranged around each 
of the dose-collecting points on the transverse plane with 
radial distance limitation |r*— p\ < ay/2 at, where the 
Gaussian tail-cutoff parameter a is normally set to 3 and 
a normalization factor is multiplied to the kernel to com- 
pensate the ignored tail contributions. 2 

The PB algorithm accommodates the density hetero- 
geneity by involving the ray-tracing integral Eq. ^ to 
derive the terma and the kernel within the convolution 
integral Eq. ([3]). The multiple integration will, however, 
increase the computational amount severely. 



C. The grid-dose-spreading algorithm 

In treatment planning, the dose grids must be fine 
enough to show dose variation in the patient with grid 
spacing as small as at or less and should be also able to 
represent distributions of any quantities. In the grid- 
dose-spreading (GDS) algorithm, the termas and the 
spreads in Eqs. (|U) and (0 are calculated at all the dose 
grids and stored in three-dimensional arrays, 



T, 



T(n) 



(6) 
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for grid i located at fi — (xi,yi,Zi) in the grid-based 
coordinate system. The number of the ray-tracing inte- 
grals is minimized by extracting out of the convolution 
integral. 



The gridded distributions, however, are not directly 
applicable to the convolution due to the coplanar con- 
straint between terma emission and dose collection in 
the PB model because the PB axis is generally angled 
to the x, y, and z grid axes with direction cosine vector 
v = (v x ,v y ,v z ). In order to resolve this difficulty, the 
disk-shaped kernel in Eq. is deformed to the best ap- 
proximate ellipsoidal kernel of the product of three Gaus- 
sian functions, 
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where q — r — p — {q x ,q y ,q z ) is the displacement vector 
from the terma-emitting to dose-collecting points and cr x , 
(T y , and a z are the grid-axial projections of the spread, 

4iP) = J qlh(p,p + q)d 3 q 

= (l-v 2 k (p)) a 2 t {p) (fc = x,y,z), (9) 

derived with the planar point kernel in Eq. ([5]). Apply- 
ing the gridded distributions, the convolution in the PB 
algorithm in Eq. (J3]) is rewritten to 



D 3 = 
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hki(qkij) 



(10) 



where Dj is the dose collected at grid j , q\j is the displace- 
ment between grids i and j, grid-axial spreading function 
hki{qkij) is the dose fraction transfered into width 5k at 
grid j, and dose-collection acceptance tui compensates 
the ignored Gaussian tails by cutoff parameter a for the 
summation. The grid-axial spreading functions and the 
acceptances are analytically given by 



hki(q) 




(11) 



(12) 



which can be quickly computed with the standard math 
library. 

Since the effective kernel volume is conserved in the de- 
formation, the computational amount is proportional to 
the cutoff cross section of the pencil beam or roughly to 
a 2 . The computational efficiency is maximized by adopt- 
ing the convolution scheme so-called "the interaction 
point of view"fi£ In this case, the acceptance-corrected 
terma Ti/(e xi e Yi e zi ) is calculated at each emitter grid i 
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and then is distributed with the fractions J} fc /^((fty)' 8 
to the nearby grids j's. The gridded dose distribution Dj 
will be formed when all the terma emissions have been 
processed as above. 

It is noted that the GDS and PB algorithms may share 
all the physical and computational models and their in- 
accuracies except for the grid quantization and the kernel 
deformation. In cases where these additional inaccuracies 
are substantially smaller than the other ones, the GDS 
and PB algorithms will be practically equivalent. 



D. Implementation to treatment planning system 
for carbon-ion radiotherapy 

Carbon-ion radiotherapy has been practiced at Na- 
tional Institute of Radiological Sciences since 1994 with 
accelerator complex HIM AC,— and original treatment 
planning system HIPLAN ; 18 ' 19 where the BB algorithm 
has been consistently used to avoid disturbance to the 
ongoing clinical studies even though Petti and Kohno 
et al. explicitly showed that the BB algorithm involves 
principled inaccuracy due to lack of beam blurring in the 
fields 

The GDS algorithm as described in Sec. Ill Cl was imple- 
mented to HIPLAN using the existing framework of the 
BB code. The terma distribution T(r) is calculated by 
applying er t — > in Eq. (JTJ) , where the depth-dose curve 
Dsq{w) is from the beam data library of HIPLAN^ for 
the range-modulated carbon-ion beams including rela- 
tive biological effectiveness (RBE) correction Invariant 
<7t = 4 mm and a = 3 are used to preserve the penumbra 
behavior of the BB algorithm for the 400-MeV/nucleon 
beams with the multileaf collimator. 

Generally, two algorithms can be impartially compared 
only under exactly the same condition except for the es- 
sential algorithmic differences in implementation. In this 
regard, we can accurately compare the GDS and BB al- 
gorithms by applying them to the identical plan on the 
single HIPLAN system. 



E. Analytic proton pencil-beam model 



Bragg curve in water, Dbb(w) in our notation, as 
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in MeVg _1 cm 2 per incident proton^ where T> n (u) is 
the parabolic cylinder function^ and Rg and <tr are 
the beam range and the range straggling in g/cm 2 re- 
lated to incident energy Eq in MeV with formulas Rq = 
0.0022 Eq 1 - 77 and a R = 0.012 R - 935 , respectively^ 

Hong et al. tabulated the projected transverse scatter 
of protons in water j 2 - which is approximated by function 
Uo{w) = 0.023 w (0.83 w/Ro + 0.17) with accuracy better 
than a fraction of a millimeter. For a pencil beam gen- 
erated at r*o with projected angular spread 9q, the diver- 
gent term 0q s proportional to distance s is quadratically 
added to the scatter, 

a t 2 ( S ) = 0.023 2 w 2 (s) fo.83 ^ + 0.17^ + 9 2 s\ 

(14) 

where w(s) — J p(tq + s' v) ds' is the water-equivalent 
depth from tq. The dose distribution for a proton pencil 
beam in water, D(r), is then given by 



D{r) = D BB (w(r}) 
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where s(r) = v ■ (r — fo) is the projected distance from 
r*o to r along v. 

In applying the GDS algorithm to this system, the 
gridded spread is derived by <j tj = a t (s(fi)), while the 
gridded terma Ti is calculated according to the definition 
as total energy released per mass in the grid-z voxel, 
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Dbb (w(s)) ds, 



(16) 



where Vq is the volume of the voxel, and S{ ni and s ou ti are 
the distances on the beam axis to enter and to exit from 
the voxel, respectively, since Eq. ((4J is not applicable to 
the infinitesimal field. The subsequent formulation in 
Sec. Ill CI is then applicable to this system. 



Since the GDS algorithm is a variant of the PB algo- 
rithms with additional approximations, it is necessary 
and sufficient for the accuracy evaluation to examine 
how the GDS calculation reproduces a modeled pencil 
beam under realistic conditions in grid spacing, trans- 
verse spread, and incident angle. For the PB model, a 
proton beam in water, or the proton pencil kernel itself, 
may be the most appropriate because the spread may be 
the largest among the ion species and the errors will be 
the clearest for simplicity. 

Bortfeld established an analytic model for the proton 



III. RESULTS 

A. Performance in carbon-ion radiotherapy 

Figure [1] shows the clinical, or RBE- weighted, dose dis- 
tributions of the GDS and BB calculations on HIPLAN 
for a clinical case of carbon-ion radiotherapy for prostate. 
The clinical target volume (CTV) consisting of the 
prostate and the seminal vesicles, the rectum as an or- 
gan at risk, and the planning target volume (PTV) with 
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FIG. 1: Clinical dose distributions in a transaxial plane (43 
cm x 28 cm) from a carbon-ion beam for prostate treatment 
calculated with the (a) BB and (b) GDS algorithms, where a 
prostate and seminal vesicles (gray lines) are included in the 
target (light gray) that partly overlaps with a rectum (dark 
gray) in a patient (gray), a horizontal beam is incident from 
the patient's left (the figure's right), and the doses are in 
linear gray scale (black for zero to white for the maximum). 

5-mm margin to the rectum side and 10-mm margin else- 
where added to the CTV, were manually segmented; 24 
The effective density distribution for heterogeneity cor- 
rection was derived from the planning CT image, 25 
with grid spacings of 1.758 mm along the right-left and 
anterior-posterior axes and 2.500 mm along the inferior- 
superior axis, which are shared by the dose distributions. 

A horizontal beam was conformed to the PTV with 
minimum 6 mm field margin using a multileaf collimator 
and was customized with a range compensator, a sculp- 
tured plastic object attached to the port, to absorb extra 
penetration of the carbon ions beyond the PTV with 3 
mm depth margin. The compensator was designed in the 
3x3 mm 2 -sized pixel-array format with steepness lim- 
ited by the maximum depth step of 15 mm considering 
the tapered structure of the milling tool. In this example, 
the rectum side of the PTV is almost parallel to the hori- 
zontal beam and the range compensation results in steep 
variation in beam range or so-called range discontinuity. 

The BB calculation in Fig. [lja) exhibits unphysically 
too sharp dose gradient at the range discontinuity around 
the rectum at risk, which could be influential on the clin- 
ical plan review. The ripples and the spikes of the range 
surface came from incomplete range compensation within 
the 3 x 3- mm 2 pixels. These artifacts have been naturally 
smeared out in the GDS calculation in Fig. [ljb). 

While the BB algorithm was designed to reproduce 
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(a) Analytic Model (b) Grid Dose Spreading 

FIG. 2: Dose distributions from a proton pencil beam in water 
projected onto the x-z plane, (a) the analytic beam model 
and (b) the corresponding GDS calculation. The z axis is 
the vertical height from the water level and the x axis is the 
relative horizontal position. The 20% and 50% isodose lines 
relative to the analytic maximum of 29.2 MeVg - cm per 
incident proton are drawn with the gray scale images. The 
embedded images show the point-spread functions. 



the error function of at = 4 mm in the penumbra re- 
gion, there was submillimeter-level disagreement between 
the BB and GDS calculations in field edge defined by 
50%-dose position, which is consistent with the grid- 
quantization error of the GDS algorithm. The BB and 
GDS calculations for the prostate treatment case took 48 
s and 66 s with SGI® Octane workstation, respectively. 
Namely, the GDS calculation was 1.4 times slower than 
the BB calculation in this case. 



B. Accuracy for angled proton pencil beam 

Figure [2] shows dose distributions in water projected 
onto the x-z plane, where a point-like 150-MeV pro- 
ton beam with projected angular spread 9 = 10 mrad 
is generated at 10 cm above water level with zenith 
angle 30°, namely with r — (lO/y 7 ^, 0, 10) cm and 
v = (— 1/2, 0, — V3/2) in the grid-based coordinate sys- 
tem with origin defined at the beam entrance point into 
water. The grid spacings in the GDS calculation are all 
2 mm along the three axes. 

At relatively shallow depth, the GDS and the ana- 
lytic model calculations are consistent within the grid 
resolution especially in the 20% isodose line while the 
50% isodose line suffers from small dose errors under the 
low-dose-gradient condition. At the Bragg peak, there is 
substantial disagreement in the 20% isodose line within 
about 3 mm. 

The embedded images in Fig. [5] show the point- 
spreading functions with projected transverse spread of 
at = 4.5 mm at the Bragg peak, where the spreading in 
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the analytic model is confined in the transverse plane and 
that in the GDS algorithm forms an ellipsoidal volume. 
In other words, the planar spreading is approximately re- 
solved into the three uncorrelated orthogonal spreading 
by cr x = 3.9 mm, a y — 4.5 mm, and a z — 2.3 mm, which 
has deformed the point-spreading function, or the point 
kernel, and consequently the dose distribution. 



IV. DISCUSSION 
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A. Speed 
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The PB algorithms superpose the dose distributions 
from the pencil beams in the field, where time-consuming 
ray-tracing integrals will severely slow down the calcula- 
tion if performed within the convolution. The algorith- 
mic implementation from the interaction point of view^ 
where multiple dose depositions are associated with each 
terma emission, may minimize the ray-tracing integrals. 
However, the dose transfer to the nearby dose grids in- 
volves relatively complicated coordinate transformations 
within the convolution, which may be still substantially 
time consuming^ 

In the GDS algorithm, the superposition is done in 
the terma space rather than in the dose space, the ray- 
tracing integrals are completely decoupled from the con- 
volution, and the geometrical emitter-collector associa- 
tions are maximally simplified to the grid-to-grid basis, 
which in all make the GDS algorithm very efficient. In 
fact, Hong et al. studied the BB and PB algorithms for 
proton radiotherapy and found that the PB calculation 
was slower than the BB calculation by factor of 73 with 
a test case<2 The corresponding slowing factor for the 
GDS algorithm was only 1.4 in our prostate treatment 
case, indicating significant speed improvement over the 
conventional PB algorithms. 

B. Accuracy 

In the GDS algorithm, the terma distribution is cal- 
culated as an intermediate quantity with grid quantiza- 
tion errors, which will propagate to the dose distribution. 
However, the quantization error will be usually negligi- 
ble with sufficiently fine grid spacing. For example, in 
the prostate treatment case with grid spacing of 1.758 
mm, the rms error from the quantization will be as small 
as 1.758/\/l2 w 0.5 mm, which is by far better than the 
realistic accuracy in target delineation. 

The kernel deformation observed in the proton pencil 
beam case will greatly depend on beam direction with 
respect to the grid axes. When the beam is angled to 
all the three grid axes, the transverse planar spreading 
is approximated by volumetric spreading that includes 
an artifactual longitudinal component. The longitudinal 
spread a[ and the transverse spread a' t of the deformed 
kernel are derived from reprojection and conservation of 



FIG. 3: The longitudinal (a[; sofid fine) and transverse (a' t ; 
dashed fine) spreads of a copfanar beam reprojected from the 
deformed kernef in the GDS afgorithm refative to the original 
transverse spread (<r t ) as a function of gantry angie (cj>). 

the spread squared, 

°[ = j£°K = yjl-v*-v$-vi a t (17) 



which are a measure of inaccuracy in distal fall off and 
a measure of accuracy in lateral penumbra, respectively, 
and will be both ^/(2/3) a t in the worst case with direc- 
tion v = (±1, ±1, ±l)/\/3. Figure [3] shows the repro- 
jected spreads as a function of gantry angle <j> m copla- 
nar beam arrangement with v y ~ 0, which explains the 
deformation of the proton pencil beam at <f> = 30° in 

sec. mm 

The spread of a pristine pencil beam limits the granu- 
larity of the dose distribution and may naturally approx- 
imate the necessary and sufficient spatial resolution for 
beam control and dose evaluation. The finer structure 
beiow the resolution and the various spatial uncertain- 
ties are normally tolerated with appropriate margins. In 
fact, the PTV should include substantial margins against 
patient setup error and internal organ motion of typi- 
cally a few to several millimeters^ for example 5 to 10 
mm for the prostate in Sec. IIII Al For clinical proton 
beams, the 20%-80% penumbra size may grow as large 
as 10 mm^ or cr t w 10/1.68 ~ 6 mm, against which, 
a field margin of 1.5 at ~ 2 at ~ 10 mm around the 
PTV is usually added. Then, even with the worst di- 
rection v — (± 1 , ± 1 , ± 1 ) / a/3, the artifact as large as 
a[ = y/{2/Z)at ~ 5 mm will be mostly covered up by 
those margins. 

Generally for broad beams, the systematic deformation 
of the kernels uniformly distributed in the field will be 
mostly compensated except for field edges in analogy to 
the kernel-tilting approximation for photon beams— The 
spread and hence the deformation will be even smaller 
with heavier-ion beams. In addition, when a vertical 
or horizontal coplanar beam is used in conjunction with 
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planning CT in treatment position, namely with <f> = or 
90 in Fig. [31 the artifact will be completely absent, which 
has been almost always the case in carbon-ion radiother- 
apy with HIMAC and will be as well with its planned 
successors^ 

There are also other sources of inaccuracy in the PB 
model^S as well as ones in real clinical systems such as 
uncertainties in patient modeling, patient setup, organ 
motion, physiological variation, and beam control. Since 
there are always requirements for practicalness in treat- 
ment planning, accuracy is important but not the abso- 
lute measure of dose calculation and should be assessed 
in a relative manner. 



C. Applicability 

The implementation of the GDS algorithm to the treat- 
ment planning system HIPLAN should have proved its 
simplicity, effectiveness, usability, and robustness neces- 
sary for dose distribution calculation in treatment plan- 
ning practice and may also indicate possibilities for 
broader applications in the future. 

For example, the framework for the single pencil beam 
in Sec. Ill El can be extended to scanned beams, where 
the pencil beams are dynamically modulated in inten- 
sity, position, range, and spread i 5 ' 7 ' 10 The instantaneous 
pencil beam, or spot beam j, is characterized by dose-per- 
fluence-per-incident-particle distribution D^Bj (w) , num- 
ber of incident particles Nj, and transverse spread <7tj,. 
at grid i. The effective terma Ti and the effective trans- 
verse spread Oti at grid i are calculated by dose- weighted 
superposition of the spot beams, 



where parameter Sj is the distance along spot beam j 
and . and s< t . are the distances to enter and to exit 
from the grid-i voxel. These gridded distributions will be 
calculated efficiently from the interaction point of view, 12 
for the GDS algorithm in Sec. Ill CI In other words, in 
the GDS algorithm, summation of the pencil beams per 
field is efficiently made in the terma and spread spaces 
and then a single volumetric convolution is also efficiently 
applied to form the dose distribution. 

Application to intensity-modulated particle-beam 
therapy to achieve uniform target dose with multidi- 
rectional nonuniform fields , 10 i 14 however, requires some 
caution because the kernel deformation of the GDS al- 
gorithm is orientation dependent and, thus, will not be 
compensated in the summed plan-dose distribution even 
in the middle of the treated volume. Appropriateness of 
the beam modulation will have to be assessed, or oth- 
erwise compromised, to be less sensitive to not only the 
kernel deformation but also all the other uncertainties 
involved in the real clinical system^ 



V. CONCLUSIONS 

A new variant of the PB algorithm, the GDS algo- 
rithm, is proposed for heavy charged particle radiother- 
apy with approximations of the gridded intermediate dis- 
tributions and a modified convolution kernel for grid-to- 
grid energy transfer. The resultant high-speed nature 
and easiness of implementation are distinctive features 
of the GDS algorithm. 

When the beam incidence is angled to all the dose-grid 
axes, the approximation will cause deformation in dose 
distribution within the size of the pristine spread. Such 
inaccuracy will have to be assessed relatively against the 
clinical tolerances and the other sources of errors in prac- 
tical situations. 
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